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SHEAR-FLEXIBLE FINITE-ELEMENT MODELS OF LAMINATED 
COMPOSITE PLATES AND SHELLS 

Ahmed K. Noor* and Michael D. Mathers* 

Langley Research Center 

SUMMARY 

Several finite-element models are applied to the linear static, stability, and vibration analy- 
sis of laminated composite plates and shells. The study is based on linear shallov^-shell theory, 
with the effects of shear deformation, anisotropic material behavior, and bending-extensional 
coupling included. Both stiffness (displacement) and mixed finite-element models are considered. 
Discussion is focused on the effects of shear deformation and anisotropic material behavior on 
the accuracy and convergence of different finite-element models. Numerical studies are pre- 
sented which show the effects of (a) increasing the order of the approximating polynomials, 

(b) adding internal degrees of freedom, and (c) using derivatives of generalized displacements 
as nodal parameters. 


INTRODUCT ION 

Although the finite-element analysis of isotropic plates and shells has received considerable 
attention in the literature, investigations of laminated composite plates and shells are rather 
limited in extent. The reliable prediction of the response characteristics of high-modulus 
fibrous composite plates and shells often requires inclusion of the transverse shear effects in 
their mathematical models. Tliis fact has been amply documented for linear static, stability, 
and dynamic problems. (See, for example, refs. 1 to 5.) 

At present there are three approaches for developing plate and shell finite-element models 
which account for shear deformation. The first approach is based on the use of three- 
dimensional isoparametric sohd elements which automatically include the shear-distortion mecha- 
nism (refs. 6 and 7). The second approach employs two-dimensional elements used with inde- 
pendent shape (or interpolation) functions for displacements and rotations (refs. 8 and 9). The 
third approach is based on the addition of effects of shear deformation to two-dimensional 
classical plate or shell elements through the use of equilibrium equations (refs. 10 and 11). 
Although it is desirable to have an element which gives accurate results regardless of how 
important the shear deformation is, most of the existing elements do not satisfy this 
requirement. 
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In the context , of the stiffness method, the first approach has the major disadvantage 
that it leads to a stiffness matrix which is (1) very large for laminated composites consisting 
of many layers and (2) highly ill conditioned for thin plates or shells. If low-order interpola- 
tion polynomials are used, the second approach leads to overly stiff elements for very thin 
plates and shells. Although the aforementioned drawbacks have been recognized and some 
improvements have been suggested, the difficulties have not been overcome. (See, e.g., 
refs. 12 to 17.) The range of validity of the third approach has not been explored. Since 
the second approach provides flexibility and simplicity in fulfilling the interelement compati- 
bihty conditions and does not result in as large a stiffness matrix as in the first approach, it 
was adopted in the present study. 

The first objective of this paper is to assess the relative merits of a number of displace- 
ment and mixed shear-flexible finite elements when applied to the linear static, stability, and 
vibration problems of laminated plates and shells. Emphasis is focused on the effects of shear 
deformation and anisotropic material behavior on the accuracy and convergence of the different 
models. The second objective is to study the effects of increasing the order of approximating 
polynomials, adding internal degrees of freedom, and using derivatives of generalized displace- 
ments as nodal parameters on the accuracy and rate of convergence of the different models. 

To the authors’ knowledge no publication exists in which the aforementioned effects are 
studied in any detail. 

The analytical formulation is based on a form of the shallow-shell theory modified to 
include the effects of shear deformation and rotary inertia. Indicial notation is used through- 
out this paper since it is particularly useful in identifying the symmetries and, consequently, 
simplifies the element development. Both triangular and quadrilateral elements are considered. 
The elements are conforming and satisfy continuity requirements of the type (continuity 
of the fundamental unknowns). 


SY MBOLS AND NOTATION 


^ajiyp’^aliyp^ 


shell compliance coefficients, inverse of shell stiffnesses 


a side length of plate or shallow shell 

^a^yp extensional stiffnesses of shell 

^a3]33 transverse shear stiffnesses of shell 

(k) (k) 

^a^yp’^a3^3 stiffness coefficients of kth layer of shell 
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^aPyp 

E 


Ef 

El,Ej 


^aPyp 

f 


Glt.Gtt 

h 

hk,hk_i 


K 


K 


IJ 


ij 


IJ 

^a/3 


'■a, a 


M 


IJ 


M, 


ap 


m 


mo,m|,m2 


portions of shell boundary over which tractions and displacements are prescribed 

bending stiffnesses of shell 

elastic modulus of isotropic materials 

error index (see eq. (36)) 

elastic moduli in direction of fibers and normal to it, respectively 
stiffness interaction coefficients of shell 
rise of shallow shells 

shear moduli in plane of fibers and normal to it, respectively 
nodal stress resultants 
local thickness of shell 

distances from reference (middle) surface to top and bottom surfaces of kth 
layer, respectively 

stiffness coefficients of shell element 

geometric or initial stress stiffness coefficients of shell element 
curvatures and twist of shell reference surface 
direction cosines, cos(x^,x^^ 
consistent mass coefficients of shell element 
bending-moment stress resultants 
number of shape functions 
density parameters of shell 
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Wi 

xtO 
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i 
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P,Pa 


Po 

Qa 

R 


s«,s'- 

IJ IJ 


T 

U 


uc 


u° 


Ush.Ua 


Uq,,W 

V 


shape or interpolation functions 

extensional (in-plane) stress resultants 

relative magnitudes of prestress components 

total number of elements in xj- or X2-direction 

total number of nodes in finite-element model 

unit outward normal to shell boundary 

consistent nodal load coefficients 

external load intensities in coordinate directions 

intensity of uniform pressure loading 

transverse shear stress resultants 

radius of curvature 

radial coordinate in circular cylindrical shell (see fig. 24 ) 

“generalized” stiffness coefficients of shell element 

kinetic energy of shell 

strain energy of shell 

complementary energy of shell 

strain energy due to prestress 

measures of shear deformation and degree of anisotropy 
displacement components in coordinate directions 
work done by internal forces 
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w,w 

i 

P 


^r>'^r+l 
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K 

X 

X 


V 

^LT 




(i) 

a. 




n,nR 


p 



<t>a 


work done by external forces 

orthogonal curvilinear coordinate system (see fig. 1) 
nodal values of 

dimensionless eigenvalues of stiffness matrix 

relative size of rth element in variable grid (eq. (37)) 

dummy coordinates of ends of rth element 

fiber orientation, angle between fiber direction and xj-axis 

constant defined in appendix D 

in-plane loading parameter 

nondimensional frequency (co\/pa^y^Ej for plates; oj\fph^~jEj for shallow 
spherical segments; coypR^jEj for circular cylinders^ 

Poisson’s ratio for isotropic materials 

Poisson’s ratio measuring strain in T-direction (transverse) due to uniaxial normal 
stress in L-direction (direction of fibers) 

natural coordinates of node i 

natural (dimensionless) coordinate system in element domain 

functionals defined in equations (1) and (2) 

density of plate or shell material 

density of kth layer of laminated shell 

uniform extensional stress* in cylindrical shell 

rotation components 


nodal displacement parameters 


a 


shell domain 


0} circular frequency of vibration of shell 


d 


a 


d 


Range of indices: 

Lowercase Latin indices 1 to m 

Uppercase Latin indices: 

I,J 1 to 5 

I,J 1 to 8 

Greek indices 1 ,2 
Finite-element-model notation: 

SQN stiffness formulation, quadrilateral element, N shape functions per fundamental 

unknown 

STN stiffness formulation, triangular element, N shape functions per fundamental 

unknown 

MQN mixed formulation, quadrilateral element, N shape functions per fundamental 

unknown 

MTN mixed formulation, triangular element, N shape functions per fundamental unknown 

SQH stiffness formulation, quadrilateral element, Hermitian interpolation functions 

MATHEMATICAL_FgRMULATI^ 

The analytical formulation is based on a form of the shallow-shell theory, with the 
effects of shear deformation, anisotropic material behavior, rotary inertia, and bending- 
extensional coupling included. (See appendix A and ref. 18.) For stability problems, the 
prebuckling stresses are assumed to be given by the momentless (membrane) theory. Two 
finite-element formulations are considered. In the first formulation (displacement model) the 
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fundamental unknowns consist of the displacement and rotation components of the shell 
reference (middle) surface, and the stiffness matrix is obtained by using Hamilton’s principle 
(which for static problems reduces to the principle of minimum potential energy). The funda- 
mental unknowns in the second formulation (mixed model) consist of the 13 shell quantities: 
generalized displacements u^j, w, and 4>ct stress resultants ^aj3> Qa- 

fig. 1 for sign convention.) The generalized stiffness matrix is obtained by using a modified 
form of the Hellinger-Reissner mixed variational principle. 

The functionals used in the development of displacement and mixed models are given by 
the following equations: 

Displacement models 

= U + U°-W-T (1) 

Mixed models 

nR(Nc|3>Ma,3>Qa-«a>w,0„) = V + U° - - W - W - T (2) 


where 


U 2 ^ J ^^o:j37pC^a^/3 ^y^p 



+ >^a|3 \<t>p'^) + ^apyp \<Pp + 

Ca3|33(^aW 3^w + 2<t>^ 3^w 


(3) 

u° 

" 2^ f ^|3'^ 



(4) 

V 

-f [Nc^(9aU/3 + ' 


Qa(<^a + aaw)]dl2 


(5) 


~ ~2 

yp ^®aj37p ^a|3 

M7P + Gq^jp ^l-yp 

^a3p3 Qa 

T = 

= (mo(uc,u^ + 

ww) + 2mjUQ,0Q, 

+ m2<^Q,0^d^2 


(7) 

W 

= y*(Po;UQ, + pw)dJ2 + 

£ H 

*'*^a 

^oi + Q/3w)n^dc 


(8) 

W 

= £ [Naisfua + Uo,) 


( + Q|3 ^-w + wjjn^dc 


(9) 


In equations (3) to (9), ^a^p> ^a|37p’ ^a^yp extensional stiffnesses, bend- 
ing stiffnesses, and stiffness interaction coefficients of the shell; transverse shear 
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stiffnesses of the shell; and \x3^3 compliance coefficients 

(see appendix A); Icq,^ are the curvature components and twist of the shell surface; XN^^g 
are the initial stress resultants (prestress field) which are proportional to the in-plane load fac- 
tor X; Pq, and p are the external load components in the orthogonal coordinate direc- 
tions Xq, and X 3 , respectively; mQ, mj, and m 2 are density parameters of the shell 
defined in appendix B; cj is the circular frequency of vibration of the shell; is the 
shell domain; and c^ are portions of the boundary over which tractions and displace- 

ments are prescribed; n^, is the unit normal to the boundary; the quantities with a tilde 
denote prescribed boundary stress resultants and displacements; and 9 q, = — . 

OXq, 

FINITE-ELEMENT DISCRETIZATION 

The shell region is decomposed into finite elements connected at appropriate 

nodes, where the superscript e refers to the element. A typical element is isolated from the 
model and the fundamental unknowns are approximated by expressions of the form; 


Displacement models 


11 

( 10 ) 

II 

( 11 ) 

4>a = 

( 12 ) 

Mixed models 


In addition to the approximations of the generalized displacements (eqs. 
stress resultants are approximated by 

( 10 ) to ( 12 )), the 


(13) 


(14) 

Qc, = 

(15) 


where superscripts identify the location and subscripts designate the ordering of nodal unknowns; 

are the shape (or interpolation) functions; ipj (i = 1 to m, J = 1 to 5) are nodal dis- 
placement parameters (including, possibly, nodeless variables); Hj (i = 1 to m, J = 1 to 8 ) 
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are nodal stress-resultant parameters; m equals the number of shape functions in the approxi- 
mation; Greek indices take the values 1,2; and a repeated lowercase Latin index denotes sum- 
mation over the range 1 to m. 

ELEMENT-BEHAVIOR REPRESENTATION 

A number of displacement and mixed finite elements having both triangular and quadri- 
lateral shapes were developed in the present study. All the elements satisfy the continuity 
conditions required by the variational principles on which they are based. Within each family 
of elements, different shape (or interpolation) functions are used for approximating the funda- 
mental unknowns. The characteristics and designations of these elements are summarized in 
table 1 and are referred to frequently in the subsequent sections. 

All the triangular elements developed are based on complete polynomial approximations 
of the fundamental unknowns, thus ensuring that the functional variation is independent of 
coordinate transformations. Most of the quadrilateral elements considered in the present study 
are of the serendipity type (refs. 19 and 20), that is, with their nodes located along the ele- 
ment boundaries. The polynomial approximations used in these elements include terms which 
are of higher order than the complete expansion, and therefore, the functional variation is 
dependent on coordinate transformation. 

In each element, the same set of shape functions is used for approximating all the fun- 
damental unknowns and the nodal parameters are selected to be the values of the fundamen- 
tal unknowns at the different nodes. However, in one of the elements (SQ8-4 element), 
polynomials of different degree were used for approximating different sets of fundamental 
unknowns (lower degree polynomials were used for approximating the rotations); in the 
SQH element, products of first-order Hermitian polynomials were chosen as shape functions 
and the nodal parameters consisted of the generalized displacements, their first derivatives, and 
mixed second derivative with respect to the dimensionless local coordinates and ^2* 

(See appendix C.) Continuity of these derivatives is enforced along the interelement bound- 
aries. Since this is not required by the variational principle, the element is overconforming. 

For the two quadrilateral stiffness elements with four and eight nodes, internal degrees 
of freedom are added through the addition of displacement modes which vanish along the 
edges of the element. Those modes are usually called bubble functions (ref. 21). The shape 
functions associated with the internal degrees of freedom are products of the equations of the 
element boundaries times another polynomial, with the product representing bubble or internal 
displacement modes (elements SQ5, SQ7, SQ9, and SQll). The case of one internal mode 
(SQ5 and SQ9) corresponds to zero degree of the latter polynomial. (See table 1 and 
appendix C.) 
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In all the elements developed, the rigid body modes that cause no straining have not 
been included explicitly in the displacement fields; rather, implicit representation of these modes 
was made. A quantitative estimate of the accuracy of rigid-body-mode representation was made 
by evaluating the six lowest eigenvalues of the element stiffness matrix. This is discussed fur- 
ther in connection with the numerical studies. 

For modeling shells with curved boundaries, isoparametric elements were used in which 
the element boundary curves are approximated by the same shape functions used in approxi- 
mating the behavior functions, that is, 

x„ = ( 16 ) 

where xj^ are the nodal values of x^,. Numerical results obtained with the use of isopara- 
metric SQ12 elements are presented in the next main section. 

FINITE-ELEMENT EQUATIONS 

The governing equations for each element are obtained by first replacing the fundamental 
unknowns by their expressions in terms of the shape functions (eqs. (10) to (15)) in the 
appropriate functional (action integral for displacement models and Hellinger-Reissner functional 
for mixed models) and then applying the stationary conditions of that functional. This leads 
to a set of equations for each element of the following form: 


Displacement models 


kJj + XK« 


if 


} — T>i 


Pj + 


(17) 


Mixed models 


and 


= 0 

IJ J IJ J 


I 

) 


(18) 


‘ ‘■i ^ 1 


where and are stiffness and geometric, or initial stress, stiffness coefficients; 


M 


ij 


IJ 


IJ 


are consistent mass coefficients; and S- are “generalized” stiffness coefficients; 


ij 


IJ 
and 


IJ 


IJ 


Pj are consistent load coefficients. The formulas for the aforementioned stiffness, mass. 
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and load coefficients are given in appendix D. For stress-analysis problems, X = co = 0; for 
bifurcation-buckling problems, co = pj = 0; and for free-vibration problems, X = Pj = 0. 

In equations (17) and (18) the range of the lowercase Latin superscripts is 1 to m; the 
range of the uppercase Latin subscripts (I,J) and (I,J) is 1 to 5 and 1 to 8, respectively. 

The K, M, and S terms are completely symmetric under the interchange of one pair of 
indices for another, each pair of indices consisting of a superscript and a subscript just beneath 

it. 

To write equations (17) and (18) in matrix form, the first superscript-subscript pair of 
each of the K, S, and M terms defines the row number and the second pair defines 
the column number. For example, in equations (17) the term k|j is located in the 
[5(i-l) + ijth row and the [5(j-l) + J]th column of the element stiffness matrix. 

In the stress-analysis problems, the internal degrees of freedom (nodal parameters associ- 
ated with bubble modes) can be eliminated without any loss of accuracy by using the static 
condensation procedure (ref. 22). In stability and vibration problems, this is not done since 
it results in approximate elemental matrices. 

The integrals in the expressions for the stiffness, mass, and load coefficients (appendix D) 
are evaluated by means of the numerical quadrature formulas presented in references 20 and 23. 
In each case, the quadrature formula selected had the least number of points required to ensure 
exact evaluation of the integrals (depending on the degree of the interpolation polynomials). 
Exceptions to this are the cases of general quadrilateral or isoparametric elements based on the 
displacement models in which the stiffness and geometric stiffness coefficients contain fractional 
rational functions that are approximated by polynomials in the numerical quadrature process. 
Each entry in the elemental matrices S and S of the mixed models (eqs. (18)) contains 
just a single term. (See appendix D.) In contrast, the entries of the matrix K of the 
displacement models (eqs. (17)) are linear combinations of aj least four terms, as implied by 
the repeated (dummy) subscripts of the coefficients K in appendix D. In view of this, the 
formation of the elemental matrices for the mixed models is simpler and was found to be less 
time consuming than for the displacement models. 

BOUNDARY CONDITIONS 

In the displacement models, only kinematic (geometric) boundary conditions need to be 
satisfied. Force (stress) boundary conditions can also be satisfied if displacement derivatives 
are chosen as nodal parameters (e.g., SQH element). The effect of introducing the stress 
boundary conditions on the accuracy of solutions is discussed in the examples in the section 
“Numerical Studies.” 
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In the mixed models, both kinematic and force (stress) boundary conditions must be 
satisfied. The boundary conditions used in the present study are listed in table 2. The 
numeral 1 in this table indicates that the nodal parameter is retained and 0 indicates that the 
nodal parameter is set to zero. 

For inclined (or curved) boundaries, it is convenient to use a modified set of nodal 
parameters including normal and tangential components of displacements and stress resultants 
at the boundary points, that is, u^,', (j>Q,', and Q^,' (see fig. 2), where 


!u, 


■') 


■a 


= C 


a,a 




r 


Vi 




(19) 


with 


r ^ 





. , 



/ 

/ = ^a,a' 

' \ 

(20) 


! 

Ma'iSV 



V 


Qa 



(21) 

^a,a 

= cos(xc^,x^/) 


(22) 


The element equations at that boundary point are modified accordingly. For example, equa- 
tions (17) are modified as follows: 




~ ^ ^ J' 


(23) 


where the relations between ^nd Kjj are given by 


Kv = ' m ' Kh 


(24) 


= e, 




a'3 a3 


(25) 


^a',i3'+3 ^a,j3+3 


(26) 
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(27) 


*^3'3' *^33 


*^ 3 ,a '+3 ^ a , a ' *^3,01+3 


(28) 


^a'+3,|3'+3 *^a+3,|3+3 


(29) 


with similar relations for 


ASSEMBLY AND SOLUTION OF EQUATIONS 

If the elemental matrices are assembled and the boundary conditions are incorporated, 
the resulting finite-element field equations can be represented in the following compact form: 


Displacement models 


K 


+ X K 


(i//) = (p) + 


(30) 


Mixed models 


Ijj"^ X Ik] 


Hi loi 

) ={■■)+ co^ 

>1 IPI 


0 0 


H 


0 ' M'U 

^ ' J 


(31) 


where [k] , [k] , j^j , and contain the stiffness, geometric stiffness, mass, and load 

distributions; and [^s] contain the “generalized” stiffness distributions; ^ and 

are the vectors of nodal unknowns composed of the subvectors i//j and Hj at the various 
nodes; and the superscript T denotes transposition. Note that in the mixed models 
(eqs. (31)), the stress resultants are assembled first. 

The matrices [k| and [s] are symmetric, positive definite, and can be banded; the 
matrices [m] and [^k] are banded symmetric; and the matrix (sj is sparse. 

For stress-analysis problems, that is, X = to = 0, the governing equations of the displace- 
ment models (eqs. (30)) can be solved by any of the efficient direct techniques published in 
the literature. (See, e.g., refs. 24 to 26.) On the other hand, the governing equations of the 
mixed models can best be solved by the hypermatrix Gaussian elimination scheme. (See 
ref. 27.) 
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For eigenvalue problems, it is convenient to modify the equations of the mixed models 
(eqs. (31)) by first eliminating the stress resultants and then rewriting the resulting equations 
in the following form: 

® + >^[k] 

where 


= (p| + [mJ |\i 


'A 


(32) 


K = [sj^ [s]-^ gj 


The matrix 



is positive definite. 


(33) 


EIGENVALUE EXTRACTION TECHNIQUES 

In the absence of the external load vector fp) , equations (30) and (31) define an alge- 
braic eigenvalue problem. For free-vibration problems X = {p} ~ 0, the natural frequencies 
are obtained by applying the subspace iteration technique presented in reference 28 to the 
equations of the displacement model. 

The technique is based on the use of simultaneous inverse iteration with Gram-Schmidt 
orthogonalization. The number of vectors used in the iteration process is more than the eigen- 
vectors required, but much less than the dimensions of the matrices considered. 

For the mixed models, the natural frequencies are obtained by applying the Sturm 
sequence technique with iterations to the modified equations (eqs. (32)). In this technique 
the desired roots are first isolated by Sturm sequence procedure, then the inverse iteration tech- 
nique is applied for the determination of individual roots along with their eigenvectors. (See 
ref. 29.) 

For bifurcation-buckling problems, where only the minimum buckling load parameter is 
required, it is more efficient to use the inverse-power method presented in reference 30 for 
both the displacement and mixed models. 

EVALUATION OF STRESS RESULTANTS 

In the mixed models, once the problem is solved, all the stress resultants are readily 
available. On the other hand, in the displacement models the stress resultants are obtained 
from the nodal displacement parameters by using the following relations: 
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Qa - ^a3l33 + WV 3 +/ 3 ) 


(35) 


The stress resultants obtained from equations (34) and (35) generally violate both the 
interior differential equilibrium and the stress-resultant boundary conditions and generate discon- 
tinuities at the element nodes. Therefore, in the present study the customary procedure of 
averaging contributions of contiguous elements at common nodes is followed. Such averaging 
is not needed for the SQH element. 

Other techniques have been suggested to improve the accuracy of the stress calculations. 
These include the integral stress technique (ref. 31), which is based on least-squares minimiza- 
tion of the stress error function within each element, and the conjugate stress method (ref. 32), 
which uses biorthogonal expansion to the displacement approximation. Both these approaches 
involve additional computational efforts and are not used in the present study. 

NUMERICAL STUDIES 

To assess the relative merits of the different displacement and mixed finite-element mod- 
els developed in this study (table 1), a large number of linear stress-analysis, free-vibration, and 
bifurcation-buckling problems are solved by these finite-element models. Particular emphasis is 
placed on the effects of shear deformation and anisotropic material behavior on the accuracy 
and rate of convergence of the different models. 

The numerical examples are aimed at clarifying a number of questions concerning each 
of the following effects on the accuracy and rate of convergence of finite-element solutions: 

(a) an increase in the order of approximating polynomials, (b) addition of internal degrees of 
freedom, and (c) use of derivatives of generalized displacements as nodal parameters. 

PLATE EVALUATION RESULTS 

Four sets of plate problems are solved which contain some of the characteristics typical 
of practical problems and at the same time are problems for which an essentially exact solu- 
tion can be obtained. In one of the problems, comparison is made with experimental results. 
The problems examined are 

(a) Stress, free vibration, and bifurcation buckling of laminated orthotropic square plates 

with simply supported edges 

(b) Stress analysis of orthotropic square plates with clamped edges 

(c) Stress and bifurcation-buckling analysis of square anisotropic plates with simply sup- 

ported edges 

(d) Stress analysis of cantilevered skew plates 
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All the models in table 1 are applied to problems (a) and (b). The higher order dis- 
placement and mixed elements are applied to problem (c). The higher order quadrilateral dis- 
placement models SQH and SQ12 are applied to problem (d). The results of these studies are 
discussed subsequently. 


Square Plates 

The first set of problems considered is that of the stress, free vibration, and bifurcation 
buckling of orthotropic and anisotropic square plates. Most of the results presented in this 
section are for the symmetrically laminated nine-layered graphite-epoxy plates shown in figure 3. 
For these plates two fiber orientations are analyzed: 

(a) Orthotropic plates with fiber orientation (0/90/0/90/0/90/0/90/0) 

(b) Anisotropic plates with fiber orientation (6 1-6 jd 1-6 16 j-d jd 1-6 Id), where 0 < 0 ^ 45° 

For orthotropic plates the total thickness of the 0° and 90° layers is the same, and for 
anisotropic plates the total thickness of the 6 and -6 layers is the same. Boundary conditions 
for both simply supported and clamped plates are considered. 

Simply Supported Orthotropic Plates 

The orthotropic plate problems are selected because an exact (analytic) solution can be 
obtained, and therefore, a reliable assessment of the accuracy of the different finite-element 
models can be made. The various solutions obtained are listed first and are discussed subse- 
quently. Since doubly symmetric deformations of the plate are considered, only one-quarter 
of the plate was analyzed, and the symmetric boundary conditions along the center line are 
listed in table 2. 

For stress-analysis problems, the plates were subjected to uniform loading p^. In addi- 
tion to studying the accuracy of the maximum displacements and stress resultants obtained by 
the various displacement and mixed models, an error index (a function of f) has been 

introduced to provide a quantitative measure of the relative accuracy of the stress resultants 
and displacements obtained by the different models. The error index is given by 


Ef = 




( 36 ) 


where 

f any of the stress resultants or generalized displacements 
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fj, fj exact and approximate values, respectively, of the function at the ith node 

KmaxI maximum absolute value of the exact function in the domain of interest (one-quarter 

of the plate) 

h total number of nodes in one-quarter of the plate 

The error index (eq. (36)) is essentially a weighted root-mean-square error. The smaller the 
error index Ef, the more accurate the approximate solution (obtained by the finite-element 
model) is. 

To study the effect of shear deformation on the performance of the different finite- 
element models, three values of the thickness ratio h/a of the plate were considered: 
h/a = 0.1, 0.01, and 0.001. As a quantitative measure of the shear deformation, the ratio of 
the strain energy due to transverse shears to the total strain energy was computed for the three 
plates. The results are shown in table 3. As can be seen from this table, the shear deforma- 
tion is quite important for the first plate and is negligible for the latter. Table 4 gives the 
values of the error index for each of the stress resultants and generalized displacements 

obtained by some of the stiffness and mixed finite-element models for two plate thicknesses 
(h/a = 0.1 and 0.01) and three different grids. An indication of the accuracy and rate of con- 
vergence of the solutions obtained by the different models is given in figures 4 and 5, and the 
effect of h/a on the accuracy of the different models is shown in figure 6 . 

The doubly symmetric free-vibration modes of the plate are analyzed by the various ele- 
ment models. An indication of the accuracy and rate of convergence of the fundamental fre- 
quency obtained by different displacement and mixed models is given in table 5 and figure 7 
for plates with thickness ratios h/a of 0.1 and 0.01. Figure 8 shows the effect of addition 
of internal degrees of freedom on the accuracy and rate of convergence of the four- and 
eight-node stiffness quadrilateral elements. Table 6 shows the rate of convergence of the three 
vibration frequencies coj 3 , C 03 and <^33 obtained by different stiffness models. 

To study the effect of the bending-extensional coupling on the accuracy of the higher 
order models, the SQ12 and SQH elements were applied to the free-vibration problem of 
two-layered orthotropic plates. Results obtained by these two elements for the two plates 
with h/a = 0.1 and 0.01 are shown in table 7 along with the exact solutions. 

As a quantitative measure of the shear deformation, the exact frequencies obtained by 
the shear-deformation and classical theories are compared in tables 5 , 6 , and 7 . 

Since the accuracy of the different elements for buckling problems is expected to be 
similar to that for vibration problems, only the SQ12 and SQH elements were applied to the 
bifurcation buckling of a plate subjected to uniaxial edge compression The results 
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obtained using a 2 X 2 grid in the plate quarter are given in table 8 along with the exact 
solutions for the three tliickness ratios h/a = 0.1, 0.01, and 0.001. 

An examination of the results obtained for simply supported orthotropic plates reveals 

(1) Although the convergence of the solutions obtained by all the displacement models 
is monotonic in character, the convergence of the lower order models is much slower than 
that of the higher order models. This is particularly true for stress resultants and for thinner 
plates. (See figs. 4 and 7.) 

(2) For the same total number of degrees of freedom, the higher order displacement 
models (e.g., SQ12 and SQH) lead to considerably more accurate results than the lower order 
models. Tliis is particularly true for stress resultants and for thinner plates. (See fig. 5.) 

The same phenomenon is observed for vibration frequencies. As an example of tliis, for plates 
with h/a = 0.1, the fundamental frequency obtained by the SQ12 and SQH elements and 
2X2 grid (corresponding to 99 and 108 degrees of freedom) agrees with the exact frequency 
to four significant digits. (See table 5.) In contrast, the error in the fundamental frequency 
obtained by the SQ4 element and 5X5 grid (108 degrees of freedom) is approximately 
2 percent. For higher frequencies and thinner plates, the accuracy of the SQ4 element dete- 
riorated much more rapidly than that of the higher order models. (See tables 5 and 6.) 

(3) The accuracy of the solutions obtained by the lower order displacement models 
(SQ4 element) is very sensitive to variations in the thickness ratio of the plate. For thinner 
plates, the accuracy of this element was found to be very poor. (See tables 4, 5, and 6.) 

This is because the assumed displacement functions require that the element edges remain 
straight, and the predominant bending deformation in thin plates is therefore poorly represented. 
This fact has been recognized by previous investigators and improvements have been suggested. 
(See, e.g., refs. 12, 14, 15, and 33.) However, no procedure exists to improve the accuracy 

of the element for all ranges of thickness ratio of the plate. 

(4) The SQ8-4 element, with different-order polynomial approximations for displacements 
and rotations, although considerably more accurate than the SQ4 element, is found to be less 
accurate than the SQ8 element. (See fig. 4.) For thin plates (h/a = 0.001), the performance 
of the SQ8-4 element was found to be unsatisfactory. (See fig. 6.) 

(5) Of all the finite-element models considered, the most accurate results for a given total 
number of degrees of freedom were obtained with the SQH element. (See fig. 5 and tables 5 
and 6.) The SQH element has the added advantage that the stress resultants are continuous 
along the interelement boundaries and no averaging is needed in their evaluation. However, in 
the presence of concentrated loads or discontinuities in the geometric or material characteristics, 
some of the nodal parameters are discontinuous and a special treatment is needed. (See, e.g., 
ref. 34.) 
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(6) Bending-extensional coupling does not appear to have any adverse effect on the accu- 
racy of the higher order displacement models. (See table 7.) 

(7) The addition of internal degrees of freedom (bubble modes) to the displacement mod- 
els results, in general, in improving the performance of the element. (See tables 4, 5, and 6 
and fig. 8.) In stress-analysis problems where the internal degrees of freedom can be eliminated 
by static condensation techniques, this is an effective way of improving the accuracy of the ele- 
ment, without affecting the accuracy of the solution. For free-vibration problems, the addition 
of internal degrees of freedom is less effective than the addition of nodes to the element. An 
exception to this is the case of the SQ8 element when applied to the analysis of higher vibra- 
tion modes of plates. In this case addition of higher order polynomial terms associated with 
internal degrees of freedom has a more pronounced effect on the accuracy than the addition 

of nodes. (Compare the frequencies obtained by SQ9 and SQ12 elements for the case m = 3, 
n = 3 in table 6.) 

(8) Whereas for the SQ4 element addition of a single internal degree of freedom results 
in considerable improvement in accuracy, for the SQ8 element three internal degrees of freedom 
have to be added before a pronounced effect on accuracy can be observed. (See fig. 8.) An 
exception to tliis is the case of higher vibration modes, where the addition of a single internal 
degree of freedom improves the accuracy of the SQ8 element substantially. (See table 6.) 

(9) The solutions obtained by the mixed models are more accurate and less sensitive to 
variations in the thickness ratio of the plate than those obtained by the displacement models 
based on the same shape functions. (See tables 4 and 5 and figs. 4, 5, and 6.) However, the 
convergence of the solutions obtained by the lower order mixed models (MT3 and MQ4) is 
slow and oscillatory in character. Also, for a given number of degrees of freedom, the accu- 
racy of the solutions obtained by mixed models is lower than that obtained by higher order 
displacement models (SQH, STIO, and SQ12). (See fig. 5.) 

Two other conclusions were found but the solutions on which they are based are not 
reported herein. These are 

(10) The accuracy of the solutions obtained by the triangular elements was found to be 
sensitive to the choice of their orientation. The best accuracy was obtained when the displace- 
ment models (ST6 and STIO) had opposite orientation to that of the mixed models (MT3 

and MT6). (See fig. 4.) The results shown in tables 4, 5, and 6 and in figures 4, 5, 6, 
and 7 were obtained for the aforementioned choice. 

(11) The effect of satisfying the force boundary conditions for the SQH element (in addi- 
tion to the kinematic conditions), was found to be insignificant. Differences occurred only in 
the fourth significant digit. 
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Before closing this section, a comparison of the elements developed in the present study 
with those previously reported in the literature is in order. Since most of the latter elements 
do not include shear deformation, the problem of an isotropic square plate with h/a = 0.01, 
for which the shear deformation is negligible, was selected. The plate had simply supported 
edges and was subjected to uniform loading p^. The convergence of solutions obtained by 
several classical plate elements was reported in reference 11. Figure 9(a), which is reproduced 
from reference 1 1 , is contrasted with figures 9(b), (c), and (d), which show the convergence 
of the center displacement w, center bending moment Mjj, and strain energy U obtained 
by a number of displacement and mixed shear-flexible elements. Except for very coarse grids 
(2 X 2 or less in the plate quarter), the higher order elements developed in the present study 
are competitive with the refined elements previously reported in the literature. The problem 
of the thin isotropic plate represents a rather severe test for the accuracy of the shear-flexible 
elements, since the accuracy of such elements reduces with the diminishing of shear deformation. 

Clamped Plates 

To study the effect of clamped edges as boundary conditions on the accuracy of the 
different stiffness models, the edges of the orthotropic plates considered in the previous sub- 
section were assumed to be totally clamped and the plates were analyzed by the different 
stiffness and mixed models. The plates were subjected to uniform loading of intensity p^. 

The standard of comparison was taken to be the solution obtained by the SQH element and 
a 6 X 6 grid in the plate quarter for h/a = 0.1, and an 8 X 8 grid for h/a = 0.01 
and 0.001. An indication of the accuracy and rate of convergence of displacements and stress 
resultants obtained by the different models is given in figure 10 for three plate thicknesses, 
namely, h/a = 0.1, 0.01, and 0.001. Also, figure 11 shows the distribution of the transverse 
displacement w and the bending moment Mjj for the thinner plates (with h/a = 0.01 
and 0.001) obtained by the higher order displacement models SQ12 and SQH and the mixed 
model MQ8 with a 2 X 2 grid in the plate quarter. As can be seen from figure 10, the 
solutions obtained by the different displacement and mixed models were, in general, less accu- 
rate than those for simply supported edges (fig. 6). This is particularly true for thinner plates. 
An exception to tliis is the SQH element, which exhibited very high accuracy and fast conver- 
gence for all thickness ratios. Also, the remarks made in the previous subsection regarding the 
effect of h/a on the accuracy and convergence of the solutions obtained by different models 
were found to apply in this case, as well. 

Anisotropic Plates 

To study the effect of anisotropy on the performance of the liigher order displacement 
models, the fiber orientations of the graphite-epoxy plate shown in figure 3 were chosen to 
be (d 1-6 Id 1-6 Id 1-6 Id l-d Id) with 0 < 6 ^ 45°. The plate had simply supported edges and 
was subjected to uniform loading of intensity p^. 
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Before the numerical studies were conducted, the effects of variations of 6 on the 
response of the plate were studied. Also, an attempt was made to introduce a quantitative 
measure of the degree of anisotropy of the plate. Since the elastic coefficients 
(with a ¥= 13) and C^pyp, ^apyp^ ^aPyp either oc = (3 and y ^ p or 

a ^ P and 7 = p) vanish for orthotropic (and isotropic) plates and are nonzero only for 
anisotropic plates, it seems reasonable to take their contribution to the total strain energy of 
the plate as a quantitative estimate of its degree of anisotropy. Henceforth, the contributions 
of the anisotropic coefficients to the total strain energy will be referred to as U^. 

Figure 12 shows the effect of variations in d on the values of the displacement w 
and the bending-moment resultant Mjj at the center of the plate as well as on the strain 
energies U, U^, and U^j^. An examination of figure 12(c) reveals that the case d = 45^ 
leads to the highest degree of anisotropy and the maximum value of the shear deformation. 
Therefore, the anisotropic plate with 6 = 45° was adopted for the convergence studies. 

An indication of the accuracy and convergence of the higher order displacement mod- 
els STIO, SQ12, and SQH and the mixed model MQ8 is given in figure 13 for the plate thick- 
nesses h/a = 0.1, 0.01, and 0.001. The standard of comparison (converged solution) was taken 
to be the solution obtained by the SQH element and an 8 X 8 grid in the whole plate. Fig- 
ure 14 shows the distribution of the transverse displacement w and the stress resultant Mjj 
for the thinner plates (h/a = 0.01 and 0.001) obtained by the SQ12 and SQH elements with 
a 4 X 4 grid, along with the converged solutions. As in the cases of simply supported and 
clamped orthotropic plates, the fastest convergence was obtained by using the SQH elements. 

The only adverse effect of the anisotropy on the performance of the elements is in the non- 
monotonic character of the convergence of stress resultants. (See fig. 13(b).) 

As a further check on the accuracy of the SQH elements in the case of anisotropic 
plates, the bifurcation-buckling problem of the eight-layered anisotropic plate shown in fig- 
ure 15 was analyzed. The plate is subjected to combined compressive and shear edge loading. 
The same plate was analyzed in reference 35 using Galerkin’s method. The results obtained 
using three grid sizes of SQH elements (in the whole plate) are given in table 9 along with 
those of reference 35. Also, the buckling mode shapes are shown in figure 15. 

Sk ew Plates 

The next problem considered is that of the stress analysis of an isotropic skew plate 
subjected to uniform transverse loading (fig. 16). The problem was selected because it includes 
a more complex set of boundary conditions and stress patterns than the ones previously 
considered. 

For this plate and these boundary conditions, an unbounded bending moment and a stress 
singularity occur at point B. (See ref. 36.) The nature of the singularity remains unaltered 
even when the shear-deformation theory (ref, 37) is used. 
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Analytical and experimental studies of this problem were reported in reference 38. The 
analytic solution was obtained by applying the mixed Hellinger-Reissner formulation in conjunc- 
tion with direct variational methods to the classical plate theory (with shear deformation 
neglected). 

The plate was analyzed with both the SQ12 and SQH elements. An indication of the 
accuracy and convergence of solutions obtained by both elements is given in figures 16(a) 
and (b). Shown in figures 16(c) and (d) are the experimental and analytical solutions of 
reference 38 compared with the present solutions. 

An examination of figures 16(c) and (d) reveals that the solutions obtained by both the 
SQH and SQ12 elements, in addition to having fast mono tonic convergence, exhibit clearly the 
sharp gradient (singularity) of the bending-moment resultant M 22 at point B. Of the two 
finite-element solutions, the SQH solution has a faster convergence and appears to be more 
accurate. Moreover, for a 4 X 4 or finer grid, the total number of degrees of freedom in 
the SQH solution is less than those in the corresponding SQ12 solution. 

SHELL EVALUATION RESULTS 

Five sets of shell problems are solved by the displacement models developed in the pres- 
ent study. Comparison is made with exact and other approximate solutions whenever available. 
These problems are 

(a) Stress and free-vibration analysis of orthotropic shallow spherical segments 

(b) Stress analysis of anisotropic shallow spherical segments 

(c) Stress analysis of an isotropic cylindrical shell with a circular cutout 

(d) Free vibrations of an orthotropic cylindrical shell 

(e) Free vibrations of an anisotropic cylindrical shell 

All the displacement models listed in table 1 are applied to problem (a). Only the 
higher order models are applied to problem (b). The isoparametric SQ12 element is applied 
to problem (c), and the SQH element is applied to problems (d) and (e). The results of 
these studies are discussed subsequently. 

Shallow Spherical Shells 

As a first application to a shallow-shell problem, consider the stress and free-vibration 
analyses of simply supported, nine-layered, graphite-epoxy spherical segments. The geometric 
and material characteristics of the shell are shown in figure 17. As for the laminated plates 
examined in the previous subsections, shallow shells with two fiber orientations have been 
analyzed: 
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(a) Orthotropic shells with fiber orientation (0/90/0/90/0/90/0/90/0) 

(b) Anisotropic shells with fiber orientation (0/-6/6/-0/0/-0/0/-0/0), with 0 < 0 ^ 45^ 

Orthotropic Shallow Shells 

For the orthotropic shells considered, analytic solutions were obtained and used as a 
standard for comparing the different finite-element solutions. Doubly symmetric deformations 
of the shell were considered, and therefore, only one-quarter of the shell was analyzed. 

For stress-analysis problems, the shells were subjected to uniform loading p^. The 
different displacement models were used to obtain solutions for three thickness ratios of the 
shell (h/a = 0.1, 0.01, and 0.001). As a quantitative measure for the shear deformation, the 
ratios of the strain energy due to transverse shear to the total strain energy of the shell were 
computed for the three shells. Results are given in table 10, and as for ortho tropic plates, 

the shear deformation is quite important for the thickest shell and is negligible for the two 

thinner shells. 

An indication of the accuracy and rate of convergence of the solutions obtained by the 
different models is given in figure 18 for the shell with h/a = 0.1. The effect of h/a on 
the accuracy of the different finite-element solutions is shown in figure 19. The distributions 
of the transverse displacement w and the stress resultants N 22 Mjj obtained by 

the higher order elements SQ12 and SQH with a 2 X 2 grid in the shell quarter are shown 
in figure 20 along with the exact solutions for the two thinner shells (h/a = 0.01 and 0.001). 

The first four doubly symmetric vibration frequencies obtained by the different displace- 
ment models are listed in table 1 1 along with the exact frequencies for two thickness ratios 
(h/a = 0.1 and h/a = 0.01). The solutions obtained using the SQ4 element were, in general, 
far removed from the exact solutions and are not reported herein. 

The orientation of the ST6 and STIO elements, for optimum accuracy, was found to be 
the same as that for orthotropic plate problems. (See fig. 4.) 

An examination of figures 18, 19, and 20 and table 11 reveals that the remarks made 
in connection with the orthotropic-plate problems regarding the effectiveness of the higher 
order models (STIO, SQ12, and SQH elements) and the effect of internal degrees of freedom, 
apply in this case as well. The apparent poor performance of the different models for the 
case of very thin shells (with h/a = 0.001) is due to the boundary -layer effects exhibited by 
the stress resultants (see fig. 20), hence the difficulties (and nonmonotonicity) in convergence 
observed in figure 19. The convergence of the total energy obtained by the higher order 
models was fast and monotonic, even for the very thin shell. (See fig. 19(d).) 
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Anisotropic Shallow Shells 


For anisotropic shells the fiber orientations were chosen to be (6 j-6 jd j-6 jd j-Q jd l-d Id) 
with 0 < 0 ^ 45°. The shells were subjected to uniform loading of intensity p^. The 
quantitative measures for the degree of anisotropy and amount of shear deformation introduced 
for anisotropic plates were used for the anisotropic shallow shells as well. 


Figure 21 shows the effect of variations in 6 on the values of the center displace- 


ment w and the center stress resultants N 
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and M 


11 


for two thickness ratios of the 


shell (h/a = 0.1 and 0.01). Also shown (fig. 21(d)) are the strain energies U, 
The maximum values of ^a/^ occur at different values of 0. 

contrasted with the anisotropic plates, for which the maximum values occurred at 


Ug, and 


U. 


sh- 

This is to be 
0 = 45°. 


The accuracy and convergence studies were conducted for shells with 0 = 45°. Fig- 
ure 22 gives an indication of the accuracy and convergence of the center displacement w and 
the strain energy U obtained by the higher order displacement models (STIO, SQ12, and SQH) 
for the three thickness ratios h/a = 0.1, 0.01, and 0.001. The standards of comparison (con- 
verged solutions) were taken to be the solutions obtained by the SQH elements. An 8 X 8 
grid was used for shells with h/a = 0.1 and 0.01, and a 10 X 10 grid was used for shells 
with h/a = 0.001. The distributions of the normal displacement w and the stress resul- 
tants N 22 ancl Mji obtained by the SQ12 and SQH elements with a 4 X 4 grid for the 
thinner shells (with h/a = 0.01 and 0.001) are shown in figure 23 along with the converged 
solutions. As in all the previous problems, the SQH solutions had the fastest convergence. 

The degradation of accuracy due to anisotropy for very thin shells, though not pronounced for 
higher order displacement models, can be clearly seen by comparing the results in figures 20 
and 23. 


Rigid Body Modes 

For shallow shells, the rigid body modes are trigonometric in character and therefore 
are only approximated by the polynomial shape functions used in the present study. To assess 
the accuracy of the approximation, the eigenvalues of the stiffness matrices of the various dis- 
placement models were computed for the three anisotropic shallow shells with h/a = 0.1, 0.01, 
and 0.001. The lowest six eigenvalues correspond to rigid body modes; the higher modes are 
straining modes. Table 12 summarizes the lowest seven eigenvalues, the maximum eigenvalues, 
and the traces of the stiffness matrices for the various models. In all cases the ratio 
was greater than 10^, which indicates that the rigid body modes are satisfactorily represented 
in these models. 
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Cylindrical Shells 


Isotropic Cylinder With a Circular Cutout 

Consider the stress analysis of an isotropic cylindrical shell with a circular cutout sub- 
jected to a uniform axial tensile stress at its free ends. The geometric characteristics of the 
shell and loading are shown in figure 24. The problem was selected to assess the accuracy of 
the isoparametric SQ12 elements in situations where high stress gradients and curved boundaries 
occur. The shell and loading are doubly symmetric, and therefore, only one-quarter of the 
shell was analyzed. 

An approximate analytic solution for the problem, assuming the cylinder to be of infi- 
nite length, was given in reference 39, where it was shown that for this shell, the shallow- 
shell approximation is valid. Therefore, the use of the SQ12 elements, with local element 
coordinates coinciding with global shell coordinates, is justified. \ difference-based variational 
solution was given in reference 40. Finite-element solutions using higher order triangular ele- 
ments were reported in reference 41. All the aforementioned solutions were based on the 
classical shell theory (with shear deformation neglected). Solution to a similar cylinder problem 
using a refined grid of shear-flexible quadrilateral elements was reported in reference 42. 

Four graded networks with 4 X 4, 5 X 4, 5 X 6, and 8X6 SQ12 elements were used 
to analyze the shell. (See fig. 25.) In an attempt to make a rational choice for the variation 
of the grid size in both the Xp and X2-directions, a variable grid parameter f was introduced 
(ref. 43 and fig. 26): 

?r = U («r+l - Ir) <37) 

where fj. is the relative size of the rth element, r? refers to each of the xp and X2“ 
coordinates, 77^. and are the coordinates of the ends of the element, and n is the 

number of elements in the ?7-direction. A second-degree polynomial variation of fj. was 
chosen, that is, 


= a + br + cr^ (38) 

where r is the element number 1 ^ r ^ n. The coefficients a, b, and c of the poly- 
nomial are determined by specifying the relative sizes of the first and last elements 
and and using the following tliree equations: 

n 

i-E ?r = 10 (39) 

" r=l 
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fj — a + b + c 


(40) 


= a + bn + cn^ (41) 

The characteristics of the grids used in the present study are shown in figure 26. 

The maximum stress concentrations and strain energies obtained by the four 

grids are given in table 13 along with results of previous investigators. Membrane stress dis- 
tributions obtained by the 4X4 and 8X6 grids are shown in figure 27. The high accuracy 
and rapid convergence of the solutions obtained by the isoparametric SQ12 elements are clearly 
demonstrated by this example. 


Orthotropic Cylinders 

The natural frequencies and mode shapes of orthotropic, two-layered, simply supported 
circular cylinders without axial restraint are studied. The problems are selected to assess the 
accuracy of the SQH elements when applied to laminated closed cylinders with high bending- 
extensional coupling. The geometric characteristics of the shells studied are shown in figure 28. 
Shells with fiber orientation (90/0) are analyzed. 

For these cylinders an analytic solution is obtained and is used as a basis for comparison 
of the finite-element solutions. It is found that for this shell, the shallow-shell (Donnell’s) 
theory approximation is valid. The doubly symmetric vibration modes of the cylinders are 
analyzed and the symmetric boundary conditions along three of the edges are applied. This 
eliminates the axial rigid body mode of the cylinder and allows obtaining the vibration modes 
having odd values of m (axial direction) and even values of n (circumferential direction). 
Initially a uniform grid with 2x2 SQH elements was used to model one octant of the cyl- 
inder (grid 1, fig. 29); however, this resulted in poor accuracy for the frequencies and mode 
shapes with n ^ 4. Subsequently, the 2X2 grid was modified to cover only one-eighth of 
the circumference (grid 2, fig. 29). This resulted in considerable improvement in the accuracy 
of the frequencies for n = 4. The frequencies obtained by the two grids are given in table 14 
along with the analytic solutions obtained by both the shear-deformation and classical shallow- 
shell theories. This table shows the decrease in accuracy as the element size-to-wavelength ratio 
increases in the circumferential direction, as indicated by the increase of n. Numerically, the 
error increases from less than 0.5 percent for m = 1, n = 2 to approximately 25 percent 
for m = 1, n = 4. The increased stiffness of the finite-element model due to the larger ele- 
ment size-to-wavelength ratio has caused a greater increase in the error of the finite-element 
analysis between the two modes. The present example shows that the SQH elements lead to 
very accurate frequencies provided the element size is less than half the wavelength of the 
vibration mode. 
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Anisotropic Cylinders 

As a final example, consider the free-vibration analysis of anisotropic two-layered circular 
cylinders. The shells have the same characteristics as those for the orthotropic cylinders dis- 
cussed in the preceding subsection, except for the fiber orientation, which is chosen to be 
(45M5). 

Solutions are obtained using tliree grids with 2 X 4, 4 X 8, and 6 X 12 SQH elements 
in the whole cylinder. (See fig. 30.) In order to eliminate the axial rigid body mode of the 
cylinder, uj is set equal to zero at the center of each grid. The fundamental frequency 
and associated mode shapes are shown in figure 31. The rapid convergence of the solutions 
obtained by the SQH elements is clearly demonstrated by this example. 

CONCLUDING R EMARKS 

Several shear-flexible finite-element models are applied to the linear static, stability, and 
vibration problems of plates and shells. The study is based on the shallow-shell theory with 
effects of shear deformation, anisotropic material behavior, and bending-extensional coupling 
included. Both stiffness (displacement) and mixed finite-element models are considered. All 
the elements examined are conforming, satisfactorily represent the rigid body modes, and 
exhibit uniform convergence for stress-analysis, free-vibration, and buckling problems. Primary 
attention in this study is given to the effects of shear deformation and anisotropic material 
behavior on the accuracy and convergence of different finite-element models. 

On the basis of the present study, the following conclusions seem to be justified: 

1 . Higher order displacement models (with cubic or bicubic interpolation polynomials) 
have the following advantages over lower order models: 

(a) Tlie total number of unknowns required for a prescribed level of accuracy is less 
in the higher order than in the lower order models. This is particularly true for stress 
resultants and for thinner plates (with negligible shear deformation). 

(b) The performance of the higher order models is considerably less sensitive to 
variations in the thickness ratio and shear deformation than that of the lower order 
models. 

2. The use of derivatives of displacements as nodal parameters (SQH element) has the 
obvious advantage that the stress resultants are defined directly at the nodes and no averaging 
is needed. In addition, this results in improving the performance of the element. However, 
in the presence of concentrated loads or discontinuities in the geometric or elastic characteris- 
tics of the shell, some of the parameters will be discontinuous and a special treatment is 
needed. 
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3. The addition of internal degrees of freedom (bubble modes) to displacement models 
results, in most cases, in improving the performance of the element. In stress-analysis problems 
where the internal degrees of freedom can be eliminated by static condensation techniques, 

this is an effective way of improving the accuracy of plate and shell elements without affect- 
ing the accuracy of the solution. For free-vibration (and buckling) problems, the addition of 
internal degrees of freedom is less effective than the addition of nodes to the element. An 
exception to this is the case of the eight-node quadrilateral element when applied to the 
analysis of higher vibration modes. In this case, addition of internal degrees of freedom has 
a much more pronounced effect on the accuracy than the addition of nodes. 

4. If mixed models are contrasted with displacement models, the following can be noted: 

(a) The development of mixed models involves considerably less algebra than the 
development of displacement models. 

(b) The performance of mixed models is, in general, insensitive to variations in the 
thickness ratio and shear deformation. 

(c) Use of lower order interpolation functions (linear or bilinear) leads to a medi- 
ocre type of performance. Considerable improvement in the performance is achieved by 
using quadratic shape functions. 

(d) For a given number of degrees of freedom, the higher order displacement mod- 
els (with cubic or bicubic interpolation polynomials) lead to higher accuracy than the 
mixed models with quadratic shape functions. The effective use of mixed models requires 
the development of efficient equation-handling techniques (e.g,, based on hypermatrix stor- 
age schemes). 

5. Whereas material anisotropy was shown to have an adverse effect on the performance 
of different displacement and mixed elements, the bending-extensional coupling does not seem 
to have any pronounced effect on the accuracy and convergence of these elements. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
November 10, 1975 
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APPENDIX A 


FUNDAMENTAL EQUATIONS OF SHEAR-DEFORMATION SHALLOW-SHELL T HEORY 
The fundamental equations of the shallow-shell theory are given in this appendix. 

STRAIN-DISPLACEMENT RELATIONSHIPS 
The relationships between strain and displacement are 
+ ^^u^) + w 

Xq;|3 ~ ^(3^ a) 


+ <t>a 

where are the extensional strains of the reference surface of the shell; Xa^ 

curvature changes and twist; and 26^,3 are the transverse shearing strain components. 

CONSTITUTIVE RELATIONS OF THE SHELL 

The relations between the stress resultants and strain components of the shell are 

^a(5 ~ ^a(5yp ^a(5yp ^p 


^a(5yp ^yp ^(x(iyp Xyp 

Qa = ^a3p3 ^^^3 
The inverse relations are given by 

^ q !|3 “ ^a(iyp ^yp ^apyp ^yp 
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APPENDIX A 


Xo|3 ®a/37p ^7p ^oc^yp ^7p 


- A^3^3 Qp 

The C, F, and D coefficients are shell stiffnesses and the A, B, and 
are shell compliances defined in appendix B. 


coefficients 
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ELASTIC COEFFICIENTS OF LAMINATED SHELLS 


ELASTIC STIFFNESSES OF THE LAYERS 

no (k) 

The nonzero stiffness coefficients *"a3|33 orthotropic layer of 

the shell referred to the directions of principal elasticity are given by 


(k) _ _ Jk)/-(k) 


^1111 “ El 7 A 


(k) 

‘=1122 


(k) p(k)/-(k) 
*^LT Et / A 


(k) 

‘=2222 


Ej y 


(k) _ ^(k) 

^n\2 ~ k^LT 


and 


(k) 

‘=1313 


G 


(k) 

LT 


(k) _ ^(k) 

‘=2323 “ kj-pL 

where the subscripts L and T denote the direction of fibers and the transverse direction, 
t>LT is Poisson’s ratio measuring the strain in the T-direction due to a uniaxial normal stress 
in the L-direction: 


'^TL El ~ *^LT Et 


X - 1 - J'lT ^TL 

and the superscript k refers to the kth layer. 
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The stiffness coefficients and Coi3^3 satisfy the following symmetry relationships 

^apyp ~ ^7paj3 “ ^payp ~ ^cxppy 


and 


‘^a3/33 " ^/33a3 “ ^3ap3 ~ ^cx33p 

If the coordinates Xq, are rotated, the elastic coefficients c^p^p and trans- 

form as components of fourth- and second-order tensors, respectively. The transformation law 
of these coefficients is expressed as follows: 

^a'^'y'p' ~ ^OL^yp ^a.,a! ^y,y' ^p,p' 


and 

^a'3|3'3 " ^a3|33 ^a,a' 

where ^(x'^'y'p' ^^e the stiffness-coefficients referred to the new coordinate 

system x^,' and 

^a,a' ~ 

ELASTIC COEFFICIENTS OF THE SHELL 


The equivalent elastic stiffnesses of the shell are given by 


NL (,,) p 2" 

\^a(5yp’ ^aj37p’ ^o:j37j^ ~ S / ^a^yp I^’ ^3’ ^3 I *^^3 

k-1 *^hk_i 


and 



where NL is the total number of layers of the shell and hj^ and hj^_| are the distances 
from the reference surface to the top and bottom surfaces of the kth layer, respectively. The 
elastic compliances of the shell Bq:^^P’ ^ajiyp^ ■^a3j33 obtained by inver- 

sion of the matrix of the elastic stiffnesses. (See ref. 18.) 
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The shell stiffnesses and compliance coefficients satisfy symmetry and transformation rela- 
tions similar to those of the stiffness coefficients of individual layers. 

The density parameters of the shell are given by 


NL ^hi, 

1^0. mj, m^ = ^ J [i, X 3 , X 3 ^J 

hk_i 


dxc 


no 

where Pg is the mass density of the kth layer of the shell. 
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SHAPE FUNCTIONS US ED IN PRESENT STUDY 
QUADRILATERAL ELEMENTS 

The expressions of the shape functions for the different elements developed in this study 
in terms of the quadrilateral coordinates given in this appendix. 


Bilinear Shape Functions 

The shape functions for the bilinear approximations (elements SQ4 and MQ4, see 
sketch (a)) are given by 


H 





where 



(with a = 1,2) are the quadrilateral coordinates of node j. 


Quadratic Shape Functions 


The shape functions for the quadratic approximations (elements SQ8 and MQ8, see 
sketch (b)) are given by 


Corner nodes 

wj 

Midside nodes 
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Cubic Shape Functions 

The shape functions for the cubic approximations (element SQ12, see sketch (c)) are 
given by 


Corner nodes 


Other nodes 

' U' - ^ 

w' = ^(l + £,£?’)(l 



Hermitian Shape Functions 


The Hermitian shape functions (element SQH, sketch (d)) used in the present study were 
products of the following set of first-order Hermite polynomials (sketch (e)): 


fl(?) = - 3? + 2) 

f2(r) = - ? + 1) 

fs(?) = - 3r - 2) 

f4(f) = - 1) 



APPENDIX C 


9v 9v 

If the order of the nodal parameters at each node is chosen to be v, 


and 


d^v 


9^1 9^2 

tions are given by 


where v denotes any of the fundamental unknowns, then the shape func- 


n' - fj({i ffte 
= fi+i'Si.ff'S: 


= fi(«i fe+i'{2 


- fi+l(Sl' fs+l^J2' 

where the subscripts i and 2 are functions of j as follows: 


j 

i 

2 ! 

... 1 

1 

1 

1 

5 

3 

1 

9 

3 

3 

13 

1 

3 


Shape Functions Associated With Nodeless Variables (Bubble Modes) 


Elements SQ5 and SQ9 


These elements have one bubble mode given by 


w" =(i - ti^Hi - 


/j = 5 for SQ5\ 
\j = 9 for SQ9/ 
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Elements SQ7 and SQll 

These elements have tliree bubble modes given by 

A 


= (l - «l^(l - £ 2 ^) 

N''^' = £,(l - £,^(1 - £ 2 ^) / 0 = 5 for SQ7; J = 9 for SQll) 


W' 


. £3(1 _ 5 ,^)(l _ J/) 




TRIANGULAR ELEMENTS 


The expressions of the shape functions for the different elements developed in this study 
in terms of the triangular (or area) coordinates 5 1 5 ^ 2 ’^3 given in the following 

sections. 


Linear Shape Functions 

The shape functions for the linear approximations (element MTS, sketch ( 0 ) in terms of 
triangular coordinates are given by 



Quadratic Shape Functions 

The shape functions for the quadratic approximations (elements ST6 and MT6, sketch (g)) 
in triangular coordinates are 

Corner nodes 

- l) (j = 2i - 1; i = 1 to 3 and is not summed) 
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Midside nodes 

= 4|| (j = 2i; i = 1 to 3 and is not summed; 



Cubic Shape Functions 

The shape functions for the cubic approximations (element STIO, sketch (h)) in triangular 
coordinates are given by 


Corner nodes (nodes 1, 4, 7) 

= i- 


^(3|i - l)(3^j - 2)^j (j = 3i - 2; i = 1 to 3 and is not summed) 



Boundary nodes 

~ T ^i ^i+l(^^i ~ 0 0 




3i - 1; nodes 2, 5, 8) , 

(^i = 1 to 3 and is not 


W-* = y li ^i+i(3^i+l - 1 ) 0 = 3i; nodes 3, 6, 9) 


summed; ^4 - 


Interior node (node 10) 

= 27$il2^3 


38 



APPENDIX D 


FORMULAS FOR COEFFICIENTS IN GOVERNING 
EQUATIONS FOR INDIVIDUA L ELEMENTS 

The expressions for the independent stiffness coefficients in equations (17) are given by 

^a3 “ .^{e) 

^a,l3+3 ^ J^is) ^“TlSp 

= X(e) IforSp ‘=«7 Hk + ‘^a3/33 9 (JnO 
4+3,3 = /"(e, f c<7/Jp kpp »yU‘ n‘ }fM>] dfl 

4+3J3+3 = [Da7(Jp dpW' + C,,3p3 w‘ mO d!l 

The independent nonzero geometric stiffness coefficients are given by 

^33 = nS/J SpN* 111 

The independent nonzero consistent mass coefficients are given by 
Ma/3 = "^0 dfi 

W 1/3+3 " ^(e) "^1 dfi 

W33 = ’^O W' ci ^2 
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'^i+3,|3+3 “ ^2 dO. 

where 5^^ is the Kronecker delta on a and j3. 

The expressions for the “generalized” stiffness coefficients in equations (18) are given by 

^a-h3-l ,7+p-l “ Jq(s) 

^a+/3-l,7+p+2 ^ J^(e) 

^a+j3+2,7+p+2 ~ ^(e) ^ 

Sa+6,|3+6 = ^a3P3 


sVl,7 = '^^(e) (^137^'' + 5a7W‘ 

C/3-1,3 = £ 


^al5- W-* dl2 


^a+j3+2,7+3 ^ '*' ^/; 


= /* 

'a+6,3 y^(e) 


Wi ac^/^-* dS2 




^a+6,/3+3 ■ ^a/3 J^(g) W-* dJ2 

The consistent nodal load coefficients are given by 


/i Ki^ 


pi 


pij = / , , /yJ p* 
3 7n(e) 


dn 


dJ2 
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In the above equations the contributions of the line integrals have been neglected for 
simplicity; k is a constant equal to 1 when a j3 and 1/2 when a = |3; the range of 
the lowercase Latin indices is 1 to m, where m is the number of shape functions; the 
range of the Greek indices is 1,2; and a repeated index denotes summation over the full 
range of the index. 

It should be mentioned that for elements with internal degrees of freedom (SQ5, SQ7, 
SQ9, and SQll), the indices ij in the expressions for Pj^, and P^ were assumed to have 
a range equal to the number of nodes in the element (i.e., 4 for SQ5 and SQ7 elements, 
and 8 for SQ9 and SQll elements). This means that the loading was distributed on the nodes 
of these elements and no loading was associated with internal degrees of freedom. 
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- 1 ^ 

ON 


TABLE L- CHARACTERISTICS OF SHEAR-FLEXIBLE FINITE-ELEMENT MODELS 

USED IN PRESENT STUDY 


Formulation 

Element 

shape 

Number 

of 

nodes 

Approximation 

Internal 

D.O.F.^ 

per 

Total number 
of D.O.F.^ 

Designation 




unknown 

Shell 

Plate*^ 

Name 

Symbol 

Displacement 

Quadrilateral 

4 

Bilinear 

None 

20 

12 

SQ4 

A 

(stiffness) 


4 


1 

25 

15 

SQ5 

+ 



4 


3 

35 

21 

SQ7 

% 



8 

Quadratic 

None 

40. 

24 

SQ8 

□ 



8 


1 

45 

27 

SQ9 




8 


3 

i 

33 

SQll 




12 

Cubic 

None 

1 

i 60 

36 

SQ12 

o 



4 

Product of first-order 

None 

1 80 

48 

SQH 

D 




Hermitian polynomials 


! 






8 for Uq,,w 
4 for 

Quadratic 

Bilinear 

^ None 

32 

( 

i 

SQ8-4 



Triangular 

6 

Quadratic 

None 

30 

■ 18 

ST6 

' 0 



10 

Cubic 


50 

30 

STIO 

7 

Mixed 

Quadrilateral 

! 4 

Bilinear 

None 

52 

32 

MQ4 

■ 



8 

Quadratic 

None 

104 

64 

MQ8 

♦ 


Triangular 

3 

Linear 

None 

39 

24 

MTS 

A 


_ _ 1 

Quadratic 


78 

48 

MT6 

• 


^Degrees of freedom. 

^Degenerate case of symmetrically laminated plates. 



TABLE 2.- BOUNDARY CONDITIONS USED IN PRESENT STUDY 


[o denotes suppressed degree of freedom; 1, free (unrestrained^ 


(a) Kinematic boundary conditions 


Boundary = Const 

u„ 

3-a 

W 

4>(x 

^-1 

3-(x 

Simple support 

1 

0 

0 

1 

0 

Qamped 

0 

0 

0 

0 

0 

1 Line of symmetry 

0 

1 

1 

0 

1 


(b) Force boundary conditions 


Boundary - Const 

N 

N, - 
3-a,3-a 

^12 

^a,0! 

M- o 
3-a,3-a 

Mi2 



Simple support 

0 


1 

0 

^a 

1 

1 

3 Kg 

Clamped 

1 

1 

1 

1 

1 

1 

1 

^IC 

*^a 

Free 

0 

1 

0 

0 

1 

0 

0 

1 

Line of symmetry 

1 

1 

0 

1 

1 

0 

0 

1 


= I for anisotropic shells and 0 for isotropic or orthotropic shells. 

a 


4 ^ 



00 


TABLE 2.- Concluded 


(c) Boundary conditions for SQH element along edge x^, = Const 


1 

Simple support 




Clamped 


Line of symmetry 

Displacement 

fa 

3cf 

a. f 

3 -a 



f 1 

a„f 

. 

; 3.5 f 

1 3-a 

■ 

f 

3^f 

^3-«f , 

^i¥ 


1 


1 



0 i 

1 

0 

1 ; 

0 

1 


1 


0 

1 

0 


1 

0 ; 

1 

0 

1 ; 

1 

0 

1 

0 

3-q: 

w 

0 

I 

1 

0 


1 

! 0 ' 

1 

1 

0 

1 ' 

1 

0 

1 

0 

1 

i 1 


, 1 

1 


, 0 : 

1 

0 

■ 1 ' 

0 

1 

0 

' 1 

1 3-q; 

i 0 

' 1 

1 

0 

1 

1 

L_ij 

1 

0 

1 

1 

0 

' 1 
1 

' 0 


^ f stands for any of the generalized displacements u^,, w, 

= 0 if force boundary conditions are imposed and 1 otherwise. 



TABLE 3.- EFFECT OF THICKNESS RATIO h/a ON TOTAL 
AND TRANSVERSE SHEAR-STRAIN ENERGIES OF PLATES 

Simply supported, nine-layered, square orthotropic plate subjected to 
uniform pressure loading p^; U denotes total strain energy of 
_ plate and denotes shear-strain energy of plate 


h/a 

UEjh/po^a^ 

Ush 

X 100 

0.1 

0.1256 

26.055 

.01 

9.2980 

.3577 

.001 

926.5123 

.0036 



2X2 

3X3 

4X4 


TABLE 4.- ERROR INDEX FOR GENERALIZED DISPLACEMENTS AND STRESS 
RESULTANTS OBTAINED BY DIFFERENT STIFFNESS AND MIXED MODELS 

Simply supported, nine- layered, square orthotropic plate subjected to unifonn loading; 


2X2 

3X3 

4X4 


212.75 

76.74 

35.60 


Values of E^ X 10 
ST6 SQ8 

f = w 


SQ 


T3 




274.91 

69.89 

33.83 


4X4 


1 

74.48 

70.82 

17.94 

2.13 

1.74 

0.40 

347.48 

213.04 

16.7 

8 

21.36 

20.82 

4.13 

.37 

.27 

.05 

116.89 

77,46 

6 . 9 ‘ 

0 

8.90 

8.76 

1.57 

.90 

.07 

.03 

44.90 

30.65 



2X2 

612.41 

58 

3X3 

289.46 

27 

4X4 

162.51 

15 


40 

928.74 

69.26 

11.08 

0.90 

219.86 

210.60 

12, 

.48 

586.61 

22.53 

2.43 

.27 

78.74 

95,72 

8 

.99 

422.01 

10.01 

.83 

.11 

36.60 

36.24 














TABLE 4.- G 


(a) Concluded 

Grid size 
(quarter 
plate) 

Values of Ej' X 10^ 

SQ4 

SQ5 

SQ7 

ST6 

SQ8 


f = 

2X2 

3X3 

4X4 

623.40 

303.48 

173.43 

598.13 

293.98 

168.99 

598.21 

293.98 

168.99 

1032.73 

713.88 

543.92 

82.22 ^ 
27.20 ! 

12.09 1 

f=Q, 

2X2 

3X3 

4X4 

1333.68 

602.92 

326.58 

1399.47 

608.81 

325.75 

1399.36 

608.47 

325.62 

1724.02 

1183.26 

891.00 

89.53 

29.53 
13.50 



f = Q 2 


2 

X 

2 

1647.63 

1741.03 

1739.31 

2800.71 

110.86 ' 

3 

X 

3 

745.45 

755.56 

754.88 

1879.48 

37,57 

4 

X 

4 

403.70 

403.65 

403.45 

1394.23 

17.40 




Grid size 
(quarter 
plate) 


SQ4 


SQ5 


SQ7 


2X2 

3X3 

4 X 4 806.46 416.11 406.94 


2X2 

3X3 

4 X 4 790.37 407.96 402.06 


2X2 

3X3 

4 X 4 781.46 404.68 399.41 

2X2 

3X3 

4 X 4 840.00 450.88 ! 445.94 


TABLE 4.- Continued 


(b) -j = 0.01 


Values of X 10 

^ for - 





ST6 

SQ8 

SQ12 

SQH 

MT3 

MQ4 


f = w 

121.58 

51.11 

9.06 

1.03 

319.96 

172.69 

27.40 

39.58 

10.51 

.54 

.19 

75.00 

66.33 

5.25 

14.22 

3.37 

.08 

.05 

37.48 

25.62 


f = <^i 1 

125.18 

38.25 

9.83 

0.97 

346.35 

209.19 

17.08 

42.43 

8.37 

.76 

.78 

108.31 

73.64 

5.67 

15.60 

2.92 

.13 

.17 

43.18 

29.71 



f = <t>^ 


124.82 

38.48 

10.27 

2 21 

365.39 

227.74 

18.37 

42.45 

' 8.97 

.83 

1.51 

112.04 

78.84 

5.87 

15.60 

3.07 

.15 

.35 

44.46 

32-41 



895.04 

260.91 

26.30 

66.51 

224.72 


28.52 

558.53 

89.56 

5-83 

17.91 

75.99 

91.84 

20.42 

403.17 

36.95 

1.76 

6.50 

44.41 








TABLE 4.- Concluded 
(b) Concluded 


Values of X 10 


2X2 

3X3 

4X4 














965.52 

309.37 

37.40 

76.21 

279.10 





689.66 

1 17.78 

9.40 

22.58 

99.72 

114.38 

859.45 

462.04 

457.63 

531.80 

52.16 

3.00 

8.83 

52.92 



106803.72 3674.66 

77899.77 1514.72 


4989.61 14280.28 14225.74 60490.19 


853.06 


1812.85 

329.35 

97.98 


838.56 

227.36 

84.28 


100.53 

80.35 

49.29 


6515.33 


180688.96 

4254.03 

2366.80 

1387.85 

97.70 

27, 

128160.13 

1638.37 

432.93 

390.99 

94.43 

38. 

97935.41 

962.99 

130.24 

148.25 

59.50 

19, 


17 

68.30 

14 

78.48 








TABLE 5.- CONVERGENCE OF MINIMUM NONDIMENSIONAL FREQUENCIES X 
OBTAINED BY DIFFERENT FINITE-ELEMENT MODELS 


Simply supported, nine-layered, square orthotropic plate;' 



Grid size 
(quarter 
plate) 

Values of X f X 1 O’ ^ for - 

h ! 

Analytic 

solution 

(a) 

SQ4 

ST6 

SQ8 

SQ9 

SQll 

STIO 

SQ12 

SQH 

MT3 

MQ4 

MT6 


2X2 

1 

1,827 

1.633 

1 1.627 

1.627 

1.626 

1.625 

1.625 

1.625 

1.587 

1.619 

1.625 

1 1.62500 

; 3 X 3 

1.714 

1.627 

1.625 

1.625 

1.625 




1.620 

1.624 


(1.88913) 1 

4X4 

1.675 

1.626 







1.623 

1.625 



S X 5 

1.657 








! 





T= 0.01 


2X2 

1 

2.010 

1.932 


1.905 

1.892 

1.896 

1.886 

1.836 

1.877 

1.887 

1.88576 

3X3 


1.939 

1.899 

1.898 

! 1.891 

1.887 

1.887 


1.879 

: 1.884 


(1.88913) 

,4X4 1 

4.232 

1.910 

1.891 

1.891 

1 

1 




1.884 

1 1.885 


• 

5 X 5 

1 

3.558 

1 

1.888 






- 

1 

1 




^Numbers in parentheses refer to classical-tlieory solutions (with both shear deformation and rotary inertia neglected). 




TABLE 6.- CONVERGENCE OF NONDIMENSIONAL FREQUENCIES X 
OBTAINED BY DIFFERENT STIFFNESS MODELS 

I Simply supported, nine-layered, square orthotropic platej] 



2X2 

3X3 

4X4 


9.197 

8.710 

8.558 


9.052 

8.542 

8.487 


(a) 
Values of 


ST6 

SQ8 

SQ9 

soil 

STIO 

SQ12 

SQH 



m = 1, 

n = 3 




6.100 

6.094 

6.088 

5.971 

5.937 

5.933 

5.930 

5.972 

5.959 

5.959 

5.931 

5.923 

5.922 


5.939 

5.934 

5.934 

1 






m = 3, 

n = 1 




6.404 

j 6.393 

6.388 

6.280 

6.249 

6.245 

6.242 

6.281 ! 

' 6.269 

6.267 

6.242 

6.235 

6.235 


6.250 

6.245 

6.245 






8.694 
8.517 
8.483 


3, n = 3 

8.559 

8.484 


8.595 

8.479 


8.768 

8.488 


8.478 


‘‘Numbers in parentheses refer to classical-theory solutions (with both shear deformation and rotary inertia 
neglected). 






TABLE 6.- Concluded 


(b) ^ = 0.01 


Values of for - 





SQ9 

SQll 

1 

STIO 

SQ12 

SQH 


Analytic 

solution 

(a) 


m = 1, n = 3 


2X2 


1.324 

1.278 

1.239 

1.159 

1.070 

1.075 

1.060 

1.03849 

3X3 

1 

1 

1.150 

1.112 

1.110 

1.072 

1.043 

1.042 


(1.0500) 

4X4 

2.885 

1.087 

1.070 

1.070 







m = 3, n = 1 


i 2X2 


1.638 

1.530 

1.499 

1.390 

1.303 

1.305 

1.292 

1.27201 1 

3X3 


1.383 

1.349 

1.347 

1.303 

1.276 

1.275 


(1.29354) ' 

4X4 

2.985 

1.318 

1.303 

1.303 





1 


m = 3, n = 3 


2X2 


2.942 

3.856 

1.975 

1.825 

1.902 

2.297 

1.703 

1.67351 

3X3 


2.242 

1.926 

1.775 

1.718 

1.709 

1.782 


(1.70022) 

4X4 

4.146 

1.935 

1.734 

1.716 







^Numbers in parentheses refer to classical-theory solutions (with both shear deformation and rotary 


inertia neglected). 


TABLE 7.- ACCURACY OF VIBRATION FREQUENCIES OBTAINED 
BY SQI2 AND SQH ELEMENTS 

Simply supported, two-layered orthotropic plate; 2X2 grid; 

^m,n = 



Values of 

[\^ n (a/h)5 X 10- 

- ^ for - 

m,n 

SQ12 

SQH 

Analytic solution 
(a) 



h/a = 

0.1 


1,1 

1.0580 

1.0578 

1.0578 

(1.1244) 

1,3 

3,1 

^ 4.8447 

4.8405 

4.8305 

(6.5730) 

3,3 

7.3000 

6.8895 

6.8757 

(9.6664) 



h/a = 

0.01 


1,1 

1.1463 

1.1303 

1.1300 

(1.1308) 

1,3 

3,1 

^ 7.1596 

6.9498 

6.7319 

(6.7644) 

3,3 

15.7435 

10.4152 

10.1118 

(10.1725) 


^Numbers in parentheses refer to classical-theory 
solutions (with both shear deformation and rotary 
inertia neglected). 



TABLE 8.- ACCURACY OF BUCKLING LOAD PARAMETER X 
OBTAINED BY SQH AND SQ12 ELEMENTS 

Simply supported, nine-layered, square ortliotropic plate subjected 
to uniaxial edge compression; N°j = -1; 2 X 2 grid in one- 
quarter of plate _ 


h 

Values of Xa^; 
1 

Ejh^ for - 

a 

SQH 

SQ12 

Analytic solution 
(a) 

0.1 

27.012 

27.014 

27.0069 (36.1597) 

.01 

36.051 

36.419 

36.0365 (36.1597) 

.001 

36.177 

69.060 

36.1585 (36.1597) 


^Numbers in parentheses refer to classical-theory 
solutions (with shear deformation neglected). 
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TABLE 9.- CONVERGENCE OF BUCKLING LOAD PARAMETER X 
OBTAINED BY SQH ELEMENTS 


Simply supported, eight-layered, anisotropic plate with fiber 
orientation (90/0/-45/45/45/-45/0/90) subjected to 
combined compressive and shear edge loadings; 

^ = 0.0072; N°, =0; N^^ = 1 ; = -1 


Grid size 
(full 
plate) 

2X2 

3X3 

4X4 


Values of Xa^/Ejh^ for — 

SQH element 

19.745 
19.194 
19.046 


Galerkin s 
method (ref. 35)^ 

19.590 


^Based on classical theory (with shear deformation 


neglected). 



TABLE 10.- EFFECT OF THICKNESS RATIO h/a ON TOTAL AND TRANSVERSE 

SHEAR-STRAIN ENERGIES OF SHELLS 

Simply supported, nine-layered, orthotropic shallow spherical shells 
(R/a = 10, f/a = 0.0125) subjected to uniform pressure loading 
Pq; U denotes total strain energy and denotes shear-strain 

_ energy 


h/a 

UETh/(po2a'>) 

Ush 

X 100 

0.1 

.01 

.001 

0.1246 

5.6983 

15.6305 

25.8512 

.2353 

.0014 


60 



TABLE 11.- CONVERGENCE OF NONDIMENSIONAL FREQUENCIES „ 

OBTAINED BY DIFFERENT STIFFNESS MODELS 

* ^ 

Simply supported, nine-layered, orthotropic shallow 

spherical shell ( ” 0.0125); 



(a) 7 = 0.1 


Grid 

size 

(quarter 

shell) 



Values of 

Xm,n X 10 

for 



Analytic 

solution 

(a) 

ST6 

SQ8 

SQ9 

SQI 1 

STIO 

SQI 2 

SQH 




m 

= 1, n - 1 





2X2 

1.639 

1.633 

1.632 

1.631 

1.631 

1.631 

1.630 

1.630 

3X3 

1.632 

1.631 

1.631 






4X4 


1.630 






(1.893) 




m 

= 1, n = 3 





2X2 

6.104 

6.096 

6.091 

5.973 

5.939 

5.935 

5.932 

5.924 

3X3 

5.974 

5.961 

5.961 






4X4 


5.936 






(10.496) 




m 

= 3, n = 1 





2X2 

6.408 

6.395 

6.390 

6.282 

; 6.251 

I 

6.246 

6.244 

6.236 

3X3 

6.283 

6.271 

6.270 






4X4 


6.247 






(12.926) 




in 

= 3, n = 3 





2X2 

9.201 

9.055 

8.695 

8.560 

8.595 

8.770 

8.478 

8.467 

3X3 

8.71 1 

8.543 

8.518 






4X4 


8.487 






(17.001) 


^Numbers in parentheses refer to classical-theory solution (with shear deformation and rotary inertia 
neglected). 



TABLE 11.- Concluded 


Values of X 10"^ for 


Analytic 

solution 


2X2 

0.2538 

0.2452 

0.2443 

3X3 

.2459 

.2422 

.2422 

4X4 


.2416 



2X2 

3X3 

4X4 


, n = 3 


(0.2414) 


2 

X 

2 

3 

X 

3 

4 

X 

4 


3, n = 3 


.986 1.833 1.909 2.303 1.709 1.680 

.783 

(1.707) 


Numbers in parentheses refer to classical-theory solution (with shear deformation and rotary inertia 
neglected). 


2X2 

2.955 

3X3 

2.250 

4X4 


















TABLE 12.- EIGENVALUES j3 OF THE STIFFNESS MATRICES FOR VARIOUS DISPLACEMENT MODELS 

Nine-layered, anisotropic shallow spherical shell with fiber 
orientation (45/-45/45/-45/45/-45/45/-45/45); R/a = 10; 
f/a = 0.0125; side length of element in coordinate directions 
is 0.25a; jS = j3a“/^Ej 



Rigid body -2.831 X lO"'^ '-3.585 X -2.388 X -1.734 X 10"’'^ -6.781 X lO"*^ -1.550 X 10''^ 4.817 X lO"'^ 

-1.458 X 10"'^ 2.821 X lO’'^ -4.073 X lO"'^ -1.166 X lO"''^ 7.991 X lO''^ 1.081 X lO-’^ 1.217 X lO''^ 

p 5.417 X 10-'^ 5.748 X 4.218 X lO''^ 1.309 X IQ-'^ ■ 1.349 X lO"'"^ 6.649 X lO"'^ 2.373 X lO"'^ 

^ 1/-I 1 A t A 1C 


_ ^6 

1.132 X 10-'^ 

1.251 X 10-'"^ 

1.338 X 10-'^ 

5.165 X 10-*5 

1.681 X 10-*'* 

1.250 X 10-*'* 

5.094 X 10- *5 

1 

^7 

4.278 X 10-5 

5.890 X 10-5 

5.890 X 10-5 

5.183 X 10-5 

1.549 X 10-5 

3.969 X 10-5 

7.870 X 10-'^ 

Straining 

^max 

9.389 

6.276 

1.175 X 10* 

o 

X 

'-n 

2.047 X 10* 

1.232 X 10* 

2.679 


Trace 

2.170 X lo' 

3.010 X lo' 

4.245 X 10* 

5.162 X 10* 

6.522 X 10* 

8.614 X 10* 

1.233 X 10* 


ON 


TABLE 12.- Concluded 


Eigenvalue 

Mode ST6 


Values of j3 for - 


F= 0.01 


-6.576 

X 

10-16 

-1.195 

X 

10-15 

-3.381 

X 

10-15 

-9.607 

X 

10-16 

-2.140 

X 

10-15 

-2.429 

X 

10-15 

-5.789 

X 

10- 

, -3.674 

X 

10-16 

-3.290 

X 

10-16 

6.559 

X 

10-17 

-5.841 

X 

10-15 

-1.312 

X 

10-15 

-3.426 

X 

10-17 

-1.665 

X 

10“ 

-1.789 

X 

10-16 

2.230 

X 

10-16 

2.532 

X 

10-16 

1.327 

X 

10-16 

-7.098 

X 

10-16 

1.631 

X 

10-15 

-5.191 

X 

10- 

2.865 

X 

10-16 

3.447 

X 

10-16 

6.021 

X 

10-16 

4.776 

X 

10-16 

-8.412 

X 

10-18 

5.699 

X 

10-15 

7.140 

X 

10- 

7.590 

X 

10-16 

2.430 

X 

10-15 

1.652 

X 

10-15 

1.852 

X 

10-15 

6.582 

X 

10-16 

7.642 

X 

10-15 

2.271 

X 

10- 

1.396 

X 

10- '5 

3.844 

X 

10-15 

3.718 

X 

10-15 

6.166 

X 

10-15 

1.817 

X 

10-15 

1.080 

X 

10-14 

9.589 

X 

10- 

6.455 

X 

10-8 

6.711 

X 

10-8 

6.711 

X 

10-8 

6.651 

X 

10-8 

4.364 

X 

10-8 

4.719 

X 

10-8 

2.091 

X 

10 - 

9.389 

X 

10-' 

6.276 

X 

10-1 

1.175 



1.175 



2.047 



1.232 



2.678 

X 

10- 

2.169 



3.007 



4.241 



5.158 



6.516 



8.608 



1.232 




0.001 


Straining 


-1.125 

X 

10-15 -1.663 

X 

10-15 

-6.682 

X 

10-17 

^ 1.006 

X 

10-15 .3 035 

X 

10-16 

-3.940 

X 

10-16 

-2.079 

X 

10- 

-5.080 

X 

10-17 .3 068 

X 

10-17 

-9.103 

X 

10-19 

-5.557 

X 

10-17 .2.039 

X 

10-16 

-2.291 

X 

10-16 

-1. 714 

X 

10- 

-1.616 

X 

10-17 _i ,11 

X 

10-17 

2.524 

X 

10-17 

-9.370 

X 

10“^^ -2.289 

X 

10-17 

-3.775 

X 

10-17 

-3.932 

X 

10- 

-5.074 

X 

10-18 1.869 

X 

10-17 

' 1.317 

X 

10-15 

1.233 

X 

10-17 .,784 

X 

10-17 

3.885 

X 

10-18 

-1.529 

X 

10- 

6.254 

X 

10-17 9.108 

X 

10-17 

3.705 

X 

10-15 

1.326 

X 

10-15 2.690 

X 

10-18 

7.217 

X 

10-17 

9,669 

X 

10- 

1.757 

X 

10-15 2.191 

X 

10-15 

3.847 

X 

10-16 

4.226 

X 

10-15 1,089 

X 

10-16 

1.327 

X 

10-15 

1.986 

X 

10- 

5.995 

X 

10-1° 1.485 

X 

10-^0 

1.485 

X 

10-10 

9.349 

X 

10-" 4.437 

X 

10-1^ 

4.730 

X 

10"^^ 

7.273 

X 

10- 

9.389 

X 

10-7 6.276 

X 

10-2 

1.175 

X 

10" ‘ 

1.175 

X 

10-1 2.047 

X 

10-^ 

1.232 

X 

10-^ 

2.679 

X 

10' 

2.167 

X 

10-1 3 007 

X 

10“^ 

4.241 

X 

10-> 

5.158 

X 

10-1 65,^ 

X 

10-^ 

8.608 

X 

10- J 

1.232 

X 

10- 



TABLE 13.- COMPARISON OF SOLUTIONS OBTAINED BY ISOPARAMETRIC SQ12 


ELEMENTS WITH THOSE OF PREVIOUS INVESTIGATORS 
[circular cylindrical shell with a circular cutout loaded in tensio^ 



Number of 

Q 

O 

1 !( ^ 

4 

Grid 

degrees of 


Membrane + 

UE/lao^h^) 

X 10 


freedom 

Membrane 

bending 

(b) 


4X4 

525 

3.643 

4.268 

4.734951 

5X4 

640 

3.691 

4.252 

4.735373 

5X6 

920 

3.712 

4.257 

4.711623 

8X6 

1415 

3.666 

4.223 

4.724903 

Finite differences 

753 

3.603 

4.096 



(ref. 40)^ 






Finite elements 

367 

3.690 

4.249 

4.729269 

(ref. 41)^ 






Analytic solution 


3.658 

4.180 



(ref. 39)2 







^Based on the classical theory (with shear deformation neglected). 
'^Strain energy in one-quarter of the shell. 



TABLE 14.- ACCURACY OF VIBRATION FREQUENCIES OBTAINED 

BY SQH ELEMENTS 

Simply supported, two-layered, orthotropic circular cy linder 
with h/R = 0.05, R/Lj = 0.5; 




Values of 

^m,n 

m,n 

SQH element 

Analytic 


Grid 1^ 

Grid 2 ^ 

solution 

(C) 

1,2 

0.5512 

— 

0.5487 (0.5494) 

1,4 

.7932 

0.6396 

.6356 ( .6473) 

3,2 

1.7173 

— 

1.7121 (1.7237) 

3,4 

! 

1.4143 

1.3390 

1.3317 (1.3581) 


^Grid 1 : 2 X 2 in shell octant. 

^Grid 2: 2 X 2 (2 elements in one-eighth of 
the circumference). 

^Numbers in parentheses refer to classical-theory 
solution (with both shear deformation and rotary inertia 
neglected). 






exact 



Number of elements per side (quarter plate) 
(a) w at point A. 


Number of elements per side (quarter plate) 
(b) <^2 point B. 


Figure 4.- Convergence of stress resultants and generalized displacements with grid refinement. Simply supported, 

nine-layered, orthotropic square plate with ^ = 0.1. 








11, exact 



Number of degrees of freedom (quarter plate) 

Figure 5.- Convergence of w and obtained by different stiffness and mixed 

models with increasii^ number of degrees of freedom. Simply supported, nine 
layered, orthotropic square plate with ^ = 0.1 and 0.01. 




exact 




SQ8-4 SQ8 SQ12 SQH 



Number of elements per side (quarter plate) 

(a) Transverse displacement w at center. 

Figure 6.- Effect of h/a on convergence of bending -moment resultant and transverse displacement w 

obtained by different stiffness and mixed models. Simply supported, nine-layered, orthotropic square plate. 







exact 





Number of elements per side (quarter plate) 



Number of elements per side (quarter plate) 


(a) 1=0.1. (b) 1=0.01. 

Figure 7.- Convergence of fundamental frequency obtained by different stiffness and mixed models with grid 
refinement. Simply supported, nine-layered, orthotropic square plate. 
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(a) ^ = 0.1; four-node quadrilaterals, 

d. 

Figure 8.- Effect of internal degrees of freedom on accuracy and convergence of vibration 
frequencies. Simply supported, nine-layered, orthotropic square plate. 
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(c) 0.01; eight-node quadrilaterals. 

3 . 

Figure 8.- Concluded. 
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Number of elements per side, n (whole plate) 

(a) w at center for previously developed (b) w at center for present 

classical-theory elements (ref. 11). shear -flexible elements. 



Number of elements per side, n (whole plate) 

(c) Mjj at center for present elements. (d) Strain energy for present 

elements. 


Figure 9.- Convergence of w, and U with grid size for shear -flexible elements 

of present study and some previously developed classical-theory elements (ref. 11). 
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Number of elements per side (quarter plate) 


(a) Transverse displacement w at center. 


Figure 10.- Effect of h/a on convergence of bendir^-moment resultant and transverse displacement w 

obtained by different stiffness and mixed models. Clamped, nine-layered, orthotropic square plate. 
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(a) Transverse displacement w at center of plate. (b) Bendii^ -moment resultant Mjj at 

center of plate. 

Figure 12.- Effect of fiber orientation on response of simply supported, nine-layered, square anisotropic 

plate subjected to uniform normal pressure loading. 
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Number of elements per side (whole plate) 

(a) Transverse displacement w at center. 

Figure 13.- Effect of h/a on convergence of bending-moment resultant 

transverse displacement w, and strain energy U obtained by higher order 
stiffness and mixed models. Simply supported, nine-layered, anisotropic 
square plate (45/-45/45/-45/45/-45/45/-45/4S). 
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Figure 15.- Buckling mode shapes for laminated anisotropic plate subjected to combined compressive and 

shear edge loadings. 





Number of elements per side (whole plate) 

(a) Transverse displacement w at point A. (b) Bending-moment resultant M22 at point B. 

Figure 16 .- Convergence of maximum transverse displacement w and bendii^-moment resultant M22 obtained 
by higher order stiffness models. Cantilevered isotropic skew plate under uniform loadii^. 
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Figure 16.- Concluded. 
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Figure 17.- Characteristics of laminated graphite -epoxy shallow shells 

used in present study. 
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(a) w at center, (b) N 22 at center. 


Figure 18.- Convergence of stress resultants, transverse displacement and strain energy with grid refinement. 
Simply supported, nine-layered, orthotropic shallow spherical segment subjected to uniform loading, ^ = 0.1. 
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Number of elements per side (quarter shell) 


(a) Transverse displacement w at center. 

Figure 19.- Effect of h/a on convergence of transverse displacement, stress 
resultants, and strain energy obtained by different stiffness models. Simply 
supported, nine-layered, orthotropic shallow -spherical segment. 
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(b) Membrane stress resultant N 22 at center. 
Figure 19.- Continued. 
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Figure 20.- Distribution of transverse displacement and stress resultants along 
center lines. Simply supported, nine-layered, orthotropic spherical segments 
^ = 0.01 and 0.001, 
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(a) Transverse displacement w at center of shell. 

Figure 21.- Effect of fiber orientation on response of a simply supported 
nine-layered, anisotropic spherical segment with a square planform 
subjected to uniform loadii^. 
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(c) Maximum bending-moment resultant M 
center of shell. 


21.- Continued. 
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Figure 22.- Effect of h/a on convergence of displacement and strain energy obtained 
by different stiffness models. Simply supported, nine-layered, anisotropic spherical 
segment with fiber orientation (45/-45/45/-45/45/-4S/45/-45/45). 
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(b) Strain energy in shell. 

Figure 22.- Concluded. 












Figure 23.- Distribution of transverse displ 


II along xg = |. 
= 0.01 and 0.001. 


Simply supported, 
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Figure 26.- Variable grid parameters. 



Figure 27.- Distribution of membrane stress concentration 
aii/OQ obtained by isoparametric SQ12 elements. 
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Figure 31.- Fundamental mode shapes and frequencies for simply supported, 
two-layered, anisotropic graphite-epoxy cylinder. 
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